library(haven)
library(emmeans)
## load data collection exercise data
## merged is a a merged long dataset of baseline and daily
dm <- as.data.table(read_sav("Merged.sav"))
1 LMM Notation
Let’s consider the formula for a relatively simple LMM:
Here as before, the i indicates the ith observation for a specific unit (e.g., person but the unit could also be classrooms, doctors, etc.) and the j indicates the jth unit (in psychology usually person).
Regression coefficients, the \(b\)s with a j subscript indicate a fixed and random effect. That is, the coefficient is allowed to vary across the units, j. As before, these coefficients in practice are decomposed into a fixed and random part:
\[
b_{0j} = \gamma_{00} + u_{0j}
\]
and we estimate in our LMM the fixed effect part, \(\gamma_{00}\), and the variance / standard deviation of the random effect or the covariance matrix if there are multiple random effects, \(\mathbf{G}\):
\[
u_{0j} \sim \mathcal{N}(0, \mathbf{G})
\]
Regression coefficients without any j subscript indicate fixed only effects, effects that do not vary across units, j. These are fixed effects and get estimated directly.
Predictors / explanatory variables, the \(x\)s with an i subscript indicate that the variable varies within a unit. Note that the outcome, \(y\)must vary within units to be used in a LMM.
In this case, the notation tells us the following:
\(y_{ij}\) the outcome variable, which varies both within and between people
\(b_{0j}\) the intercept regression coefficient, which is both a fixed and random effect
\(b_1\) the regression coefficient, slope, for the first predictor, which is a fixed effect only
\(x_{1j}\) the first predictor/explanatory variable, this is a between unit variable only, as the lack of an i subscript indicates it does not vary within units. It could not have a random slope.
\(b_2\) the regression coefficient, slope, for the second predictor, which is a fixed effect only
\(x_{2ij}\) the second predictor/explanatory variable, this variable varies within individuals as shown by its i subscript. It could have a random slope, although in this model, it only has a fixed effect slope.
\(\varepsilon_{ij}\) the residuals, these vary within and between units.
The following decision tree provides some guide to when a predictor / explanatory variable can be a fixed and random effect.
Type of effect decision tree
Let’s see two examples of putting this basic model into practice.
summary(lmer(mood ~ age + stress + (stress | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + (stress | ID)
## Data: dm
##
## REML criterion at convergence: 570.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.525 -0.510 0.188 0.588 2.361
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.7759 0.881
## stress 0.0644 0.254 -0.78
## Residual 0.9388 0.969
## Number of obs: 183, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.3282 0.9345 43.1574 6.77 2.7e-08 ***
## age -0.0131 0.0380 36.7556 -0.34 0.7329
## stress -0.2762 0.0758 27.3002 -3.64 0.0011 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.951
## stress -0.346 0.067
technically note now with two random effects, we assume that the random effects, \(u_{0j}\) and \(u_{2j}\), which we collectively denote \(\mathbf{u}_{j}\) follow a multivariate normal distribution with covariance matrix \(\mathbf{G}\).
Based on the little decision chart, between unit only variables, like \(age_j\) and \(Bstress_j\)cannot be random effects. Also, while it is technically possible for something to only be a random effect without a corresponding fixed effect, its not common and not recommended as it would be equivalent to assuming that the fixed effect, the mean of the distribution, is 0, which is rarely appropriate.
2 Interactions in LMMs
Interactions in LMMs work effectively the same way that interactions in GLMs do, although there are a few nuances in options and possible interpretations. Using the notation from above, let’s consider a few different possible interactions.
2.1 Cross Level (Between and Within Unit) Interactions
First, let’s take our model with age and stress and include an interaction. Here is the model without an interaction.
## long way
summary(lmer(mood ~ age + stress + age:stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age + stress + age:stress + (1 | ID)
## Data: dm
##
## REML criterion at convergence: 576.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.374 -0.443 0.163 0.580 2.311
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.433 0.658
## Residual 0.996 0.998
## Number of obs: 183, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.6818 2.3243 176.9101 1.15 0.25
## age 0.1476 0.1014 177.7585 1.46 0.15
## stress 0.7967 0.5968 168.1832 1.33 0.18
## age:stress -0.0480 0.0263 165.4710 -1.82 0.07 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age stress
## age -0.993
## stress -0.914 0.914
## age:stress 0.904 -0.915 -0.994
## short hand in R for simple main effect + interaction
## identical, but shorter to the above
summary(lmer(mood ~ age * stress + (1 | ID), data = dm))
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mood ~ age * stress + (1 | ID)
## Data: dm
##
## REML criterion at convergence: 576.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.374 -0.443 0.163 0.580 2.311
##
## Random effects:
## Groups Name Variance Std.Dev.
## ID (Intercept) 0.433 0.658
## Residual 0.996 0.998
## Number of obs: 183, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.6818 2.3243 176.9101 1.15 0.25
## age 0.1476 0.1014 177.7585 1.46 0.15
## stress 0.7967 0.5968 168.1832 1.33 0.18
## age:stress -0.0480 0.0263 165.4710 -1.82 0.07 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age stress
## age -0.993
## stress -0.914 0.914
## age:stress 0.904 -0.915 -0.994
The relevant, new, part is the interaction term, \(b_3\), a fixed effect in this case. If we focus just on that one term, we see that the coefficient, \(b_3\) is applied to the arithmetic product of two variables, here age and stress. As it happens, one of them, age, varies only between units whereas the other, stress, varies within units. You will sometimes see this termed as “cross level” interaction between it involves a between and within varying variable.
\[
b_3 * (age_{j} * stress_{ij})
\]
As with interactions for regular GLMs, interactions in LMMs can be interpretted in different ways. The two common interpretations are easiest to see by factoring the regression equation. Here are three equal equations that highlight different ways of viewing the interaction.
In the latter two formats, it highlights how the simple effect of stress varies by age and how the simple effect of age varies by stress.
The nuance in LMMs comes in because some variables vary only between units and others within units. For example, when interpretting the interaction with respect to the simple effect of stress, we could say that the association between daily stress and mood on the same day depends on the age of a participant. Conversely, when interpretting with respect to the simple effect of age, we could say that the association of participant age and average mood depends on how stressed someone is feeling on a given day. Age varies only between people, stress varies within people, so that must be taken into account in the interpretation.
2.2 Between Unit Interactions
The same approach would work with other type of variables in LMMs. For example, here we have a model with age and female as predictors. Both vary only between units.
When interpretting the interaction with respect to the simple effect of female, we could say that the association between participant sex and average mood depends on the age of a participant. Conversely, when interpretting with respect to the simple effect of age, we could say that the association of participant age and average mood depends on participant sex.
When interpretting the interaction with respect to the simple effect of stress, we could say that the association between daily stress and mood on the same day depends on same day energy level. Conversely, when interpretting with respect to the simple effect of energy, we could say that the association of daily energy and same day mood depends on how stressed someone is feeling on a given day.
2.4 ‘Prediction’ Interpretation of Interactions
When one variable in an interaction varies within units and particularly when it is a random effect, another way that people sometimes interpret interactions is that the other variable is ‘predicting’ the random effect. This occurs most often when the moderator variable varies between units only. An example may be clearer than words.
First, we fit a model with age and stress as fixed effects predictors and a random slope as well for stress, stored in m1. Then we fit the same model but adding a fixed effect interaction between stress (within) and age (between). Calculating the difference in the variance of the random stress slope yields a sort of \(R^2\) measure of the variance in the random slope explained by the stress x age interaction.
The corresponding R code is:
## main effects only
m1 <- lmer(mood ~ age + stress + (1 + stress | ID), data = dm, REML=FALSE)
## interaction model
m2 <- lmer(mood ~ age * stress + (1 + stress | ID), data = dm, REML=FALSE)
## summary of both models, get random slope for stress variance
summary(m1)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mood ~ age + stress + (1 + stress | ID)
## Data: dm
##
## AIC BIC logLik deviance df.resid
## 574.2 596.6 -280.1 560.2 176
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.581 -0.511 0.215 0.589 2.356
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.6643 0.815
## stress 0.0587 0.242 -0.77
## Residual 0.9388 0.969
## Number of obs: 183, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 6.2889 0.9103 44.2464 6.91 1.5e-08 ***
## age -0.0115 0.0370 37.5843 -0.31 0.75696
## stress -0.2754 0.0741 27.4601 -3.72 0.00092 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age
## age -0.951
## stress -0.346 0.068
summary(m2)
## Linear mixed model fit by maximum likelihood . t-tests use Satterthwaite's
## method [lmerModLmerTest]
## Formula: mood ~ age * stress + (1 + stress | ID)
## Data: dm
##
## AIC BIC logLik deviance df.resid
## 572.8 598.5 -278.4 556.8 175
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.628 -0.375 0.201 0.596 2.365
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## ID (Intercept) 0.5051 0.711
## stress 0.0497 0.223 -0.73
## Residual 0.9361 0.968
## Number of obs: 183, groups: ID, 50
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 2.3520 2.2723 44.8968 1.04 0.306
## age 0.1620 0.0993 45.3903 1.63 0.110
## stress 0.8907 0.6317 53.3680 1.41 0.164
## age:stress -0.0517 0.0279 54.2982 -1.85 0.069 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) age stress
## age -0.993
## stress -0.929 0.928
## age:stress 0.920 -0.931 -0.994
## variance in random stress slope explained by age
(.0587 - .0497) / .0587
## [1] 0.1533
Because age and stress interact in the model (as a fixed effect), the average stress slope no longer has to be the same for everyone. It can differ depending on their age. Specifically, the predicted average (fixed effects) stress slope for the jth person is \(b_{2j} + b_3 * age_j\). Recall that we normally break up random slopes into a fixed and random part:
\[
b_{2j} = \gamma_{20} + u_{2j}
\]
Without the interaction by age, any participant level differences must go into the \(u_{2j}\) component, on which the variance/standard deviation of the slope is calculated. Now the simple stress slope for a given participant is:
\[
\gamma_{20} + b_3 * age_j + u_{2j}
\]
so now the \(u_{2j}\) component, on which the variance/standard deviation of the slope is calculated, captures deviations from the fixed part, which includes both the average stress slope and a modification based on participant age. To the extent that \(b_3\) is different from 0, this will essentially reduce some differences that otherwise go into the \(u_{2j}\) component and thus will explain some of the variance in the random slope.
This is always true, in a way, with an interaction, but outside of LMMs we do not have random effects and so when we allow slopes to differ for different groups people, we do not know what an individual person’s slope was without the interaction so have no reference point. Put differently, outside of LMMs, in regular GLMs, we always assume that \(u_{2j} = 0\) so we wouldn’t really think about ‘predicting’ it when it is fixed at 0. In GLMs we only focus on how we allow the average slope to differ by level of the moderator. In LMMs we can interpret it the same way or we can interpret the moderator as ‘predicting’ the random slope.
3 Continuous Interactions in R
Aside from the notes about some minor interpretation differences, in general interactions in LMMs are analysed, graphed, and interpretted the same way as for GLMs.
First to avoid any issues around diagnostics etc. from haven labeled type data, we will convert the variable we are going to work with to numeric. Then we fit a LMM with an interaction between stress and neuroticism, mood as the outcome and a random intercept as the only random effect.
A quick check of the model diagnostics suggests that there may be an outlier on the residuals, but it is not too extreme and otherwise the data look fairly good. The residuals do not appear to follow a normal distribution that closely, partly due to the long left tail.
Applying transformations to left skewed data is more difficult as generally transformations work on long right tails. A solution is to reverse the variable, apply the transformation and then again reverse it so that the direction is the same as it originally was. We could try a square root transformations which is milder than a log transformation. To reverse it, we subtract the variable from the sum of its minimum and maximum. Next we take its square root, then we reverse by again subtracting from the sum of its minimum and maximum, but square root transformed.
The transformation appears to have modestly helped the distribution of residuals. Its not that clear whether it was bad enough to begin with and whether the transformation improved it enough that it is worth the difficulty in interpretation (mood is now square root transformed and that must be incorporated into its interpretation). For the lecture, we will proceed with the transformation, but in practice, consider whether this is worth it or only adds difficulty in understanding without improving / changing results much. We still see an extreme value present. We can remove that, as we have discussed in depth in previous lectures, update the model, and re-run diagnostics.
After removing one extreme value, another point is classified as an extreme value. We could remove that too. Note that its important that we remove it from the dmnoev dataset, not the original dm dataset, because the row numbers from these two datasets do not match and we are only removing the additional “new” extreme value found. It is relatively rare to actually iterate through this process of identifying and removing extreme values multiple times, but I do it to show how it woudl work if you wanted to.
After removing the additional extreme value, we would update the model again, and check diagnostics again.
Finally no extreme values are identified and the distributions are approximately normal. In practice it could take additional rounds of extreme value removal or you may decide to stop at one round.
At this point, we have a “clean” model and we can look at a summary of it. Although we have used summary() a lot in the past, we’ll introduce another function to help look at lmer() model results, modelTest(). In this lecture, we will only learn and interpret part of its output, with the rest of the output from modelTest() covered later. In addition to get nicely formatted results rather than a set of datasets containing the results, we use the APAStyler() function.
APAStyler(modelTest(mnoev))
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Term Est Type
## 1: (Intercept) 1.03*** [ 0.54, 1.53] Fixed Effects
## 2: neuroticism 0.06 [-0.01, 0.12] Fixed Effects
## 3: neuroticism:stress -0.02* [-0.04, 0.00] Fixed Effects
## 4: stress 0.08 [-0.05, 0.22] Fixed Effects
## 5: sd_(Intercept)|ID 0.19 Random Effects
## 6: sigma 0.26 Random Effects
## 7: Model DF 6 Overall Model
## 8: N (Groups) ID (50) Overall Model
## 9: N (Observations) 181 Overall Model
## 10: logLik -37.53 Overall Model
## 11: AIC 87.05 Overall Model
## 12: BIC 106.24 Overall Model
## 13: Marginal R2 0.17 Overall Model
## 14: Marginal F2 0.21 Overall Model
## 15: Conditional R2 0.45 Overall Model
## 16: Conditional F2 0.83 Overall Model
## 17: neuroticism (Fixed) 0.00/-0.01, p = .110 Effect Sizes
## 18: stress (Fixed) 0.01/ 0.01, p = .231 Effect Sizes
## 19: neuroticism:stress (Fixed) 0.03/ 0.03, p = .015 Effect Sizes
The results show the regression coefficients, asterisks for p-values, and 95% confidence intervals in brackets for the fixed effects, the standard deviations of the random effects, the model degrees of freedom, which is how many parameters were estimated in the model total, and the number of people and observations. For now, we will ignore all the output under row 9, N (Observations). In the case of this model we can see the following.
A LMM was fit with 181 observations from 50 people. There was a significant neuroticism x stress interaction, b [95% CI] = -0.02 [-.04, .00], p < .05.
We can also pass multiple model results in a list together, which puts the results side by side. This is particularly helpful for comparing models with and without covariates or to evaluate whether removing extreme values changed the results substantially.
APAStyler(list(
EV = modelTest(m),
NOEV = modelTest(mnoev) ))
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Term EV NOEV
## 1: (Intercept) 1.03*** [ 0.52, 1.54] 1.03*** [ 0.54, 1.53]
## 2: neuroticism 0.06 [-0.02, 0.13] 0.06 [-0.01, 0.12]
## 3: neuroticism:stress -0.02* [-0.04, 0.00] -0.02* [-0.04, 0.00]
## 4: stress 0.08 [-0.06, 0.22] 0.08 [-0.05, 0.22]
## 5: sd_(Intercept)|ID 0.19 0.19
## 6: sigma 0.28 0.26
## 7: Model DF 6 6
## 8: N (Groups) ID (50) ID (50)
## 9: N (Observations) 183 181
## 10: logLik -46.99 -37.53
## 11: AIC 105.97 87.05
## 12: BIC 125.23 106.24
## 13: Marginal R2 0.15 0.17
## 14: Marginal F2 0.18 0.21
## 15: Conditional R2 0.41 0.45
## 16: Conditional F2 0.69 0.83
## 17: neuroticism (Fixed) 0.00/0.00, p = .128 0.00/-0.01, p = .110
## 18: stress (Fixed) 0.01/0.01, p = .266 0.01/ 0.01, p = .231
## 19: neuroticism:stress (Fixed) 0.03/0.03, p = .024 0.03/ 0.03, p = .015
## Type
## 1: Fixed Effects
## 2: Fixed Effects
## 3: Fixed Effects
## 4: Fixed Effects
## 5: Random Effects
## 6: Random Effects
## 7: Overall Model
## 8: Overall Model
## 9: Overall Model
## 10: Overall Model
## 11: Overall Model
## 12: Overall Model
## 13: Overall Model
## 14: Overall Model
## 15: Overall Model
## 16: Overall Model
## 17: Effect Sizes
## 18: Effect Sizes
## 19: Effect Sizes
These results show us that we have effectively the same results with or without those two extreme values we omitted. Indeed the only changes really observed are small changes in the 95% confidence intervals and a slight decrease in the residual standard deviation, sigma, with the extreme values removed (outside, of course of the models being based on 181 instead of 183 people). In a case like this, it will not make any real difference to the interpretation if the model with extreme values or without extreme values is used.
We can understand what the interaction means in the same way as for GLMs. The following 3D figure shows the slight warp in the regression plane between stress and neuroticism with (transformed) mood. The warping is because of the interaction.
##
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
##
## last_plot
## The following object is masked from 'package:stats':
##
## filter
## The following object is masked from 'package:graphics':
##
## layout
3.1 Plotting
Typically, 3D graphs are not used to plot interactions. Instead, a few exemplar lines are graphed showing the slope of one variable with the outcome at different values of the moderator. As with GLMs, we can use the visreg() function. Here, we’ll use neuroticism as the moderator. A common approach to picking level of the moderator is to use the Mean - 1 SD and Mean + 1 SD. To do that, we first need the mean and standard deviation of neuroticism, which we can get using egltable() after excluding duplicates by ID, since neuroticism only varies between units.
egltable("neuroticism", data = dmnoev[!duplicated(ID)])
The results show, in an easier to interpret way, what the negative interaction coefficient of \(b = -.02\) means, people with higher levels of neuroticism are more sensitive to the effects of stress. People lower in neuroticism are relatively less sensitive to the effects of stress, although in both cases, higher stress is associated with lower (transformed) mood.
Another common way of picking some exemplar values is to use the 25th and 75th percentiles. These work particularly well for very skewed distributions where the mean +/- SD could be outside the observed range of the data. Again we exclude duplicates by ID and then use the quantile() function to get the values, 6, and 9 for the 25th and 75th percentiles.
When working with models that have interactions, a common aid to interpretation is to test the simple effects / slopes from the model. For example, previously we graphed the association between stress and transformed mood at M - 1 SD and M + 1 SD on neuroticism. However, although visually both lines appeared to have a negative slope, we do not know from the graph alone whether there is a significant association betwee stress and transformed mood at both the low (M - 1 SD) and high (M + 1 SD) levels of neuroticism. To answer that, we need to test the simple slope of stress at specific values of neuroticism. Our default model does actually give us one simple slope: it is the simple slope of stress when \(neuroticism = 0\). However, as we can tell from the mean and standard deviation of neuroticism, 0 is very far outside the plausible range of values so that simple slope given to us by default from the model is not too useful. We could either center neuroticism and re-run the model, which would get us a different simple slope, or use post hoc functions to calculate simple slopes.
We will use the emtrends() function from the emmeans package to test the simple slopes. This function also works with GLMs, for your reference.
The emtrends() function take a model as its first argument, then the variable that you want to calculate a simple slope for, here stress, the argument at requires a list of specific values of the moderator, and then we tell it how we want degrees of freedom calculated (note this only applies to lmer models). We store the results in an R object, mem and then call summary() to get a summary table. The infer = TRUE argument is needed to summary() if you want p-values.
mem <- emtrends(mnoev, var = "stress",
at = list(neuroticism = c(7.2 - 2.08, 7.2 + 2.08)),
lmer.df = "satterthwaite")
summary(mem, infer=TRUE)
The relevant parts of the output, for us, are the columsn for stress.trend which are the simple slopes, the values of neuroticism which tell us at what values of neuroticism we have calculated simple slopes, the confidence intervals, lower.CL and upper.CL, 95% by default, and the p-value. From these results, we can see that when \(neuroticism = 5.12\) there is no significant association between stress and transformed mood, but there is when \(neuroticism = 9.28\).
3.3 Sample Write Up
With all of this information, we can plan out some final steps for a polished write up of the results. First, let’s get exact p-values for all our results. We can do this through options to pcontrol in APAStyler(). We also re-print the simple slopes here.
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Term EV
## 1: (Intercept) 1.03p < .001 [ 0.52, 1.54]
## 2: neuroticism 0.06p = .133 [-0.02, 0.13]
## 3: neuroticism:stress -0.02p = .026 [-0.04, 0.00]
## 4: stress 0.08p = .273 [-0.06, 0.22]
## 5: sd_(Intercept)|ID 0.19
## 6: sigma 0.28
## 7: Model DF 6
## 8: N (Groups) ID (50)
## 9: N (Observations) 183
## 10: logLik -46.99
## 11: AIC 105.97
## 12: BIC 125.23
## 13: Marginal R2 0.15
## 14: Marginal F2 0.18
## 15: Conditional R2 0.41
## 16: Conditional F2 0.69
## 17: neuroticism (Fixed) 0.00/0.00, p = .128
## 18: stress (Fixed) 0.01/0.01, p = .266
## 19: neuroticism:stress (Fixed) 0.03/0.03, p = .024
## NOEV Type
## 1: 1.03p < .001 [ 0.54, 1.53] Fixed Effects
## 2: 0.06p = .115 [-0.01, 0.12] Fixed Effects
## 3: -0.02p = .017 [-0.04, 0.00] Fixed Effects
## 4: 0.08p = .238 [-0.05, 0.22] Fixed Effects
## 5: 0.19 Random Effects
## 6: 0.26 Random Effects
## 7: 6 Overall Model
## 8: ID (50) Overall Model
## 9: 181 Overall Model
## 10: -37.53 Overall Model
## 11: 87.05 Overall Model
## 12: 106.24 Overall Model
## 13: 0.17 Overall Model
## 14: 0.21 Overall Model
## 15: 0.45 Overall Model
## 16: 0.83 Overall Model
## 17: 0.00/-0.01, p = .110 Effect Sizes
## 18: 0.01/ 0.01, p = .231 Effect Sizes
## 19: 0.03/ 0.03, p = .015 Effect Sizes
Now we will make a polished, finalized figure. I have customized the colours, and turned off the legends. In place of legends, I have manually added text annotations including the simple slopes and confidence intervals and p-values for the simple slopes1 For your reference, it took about 6 trial and errors of different x and y values and angles to get the text to line up about right. I did not magically get the values to use to get a graph that I thought looked nice. That is why I think sometimes it is easier to add this sort of text after the fact in your slides of papers rather than building it into the code..
visreg(m, xvar = "stress",
by = "neuroticism", overlay=TRUE,
breaks = c(7.2 - 2.08, 7.20 + 2.08),
partial = FALSE, rug = FALSE, gg=TRUE,
xlab = "Daily Stress",
ylab = "Predicted Daily Mood") +
scale_color_manual(values = c("5.12" = "black", "9.28" = "grey70")) +
theme_pubr() +
guides(colour = FALSE, fill = FALSE) +
annotate(geom = "text", x = 4.2, y = 1.21, label = "High Neuroticism: b = -.03 [-.09, .02], p = .26",
angle = -10) +
annotate(geom = "text", x = 4.4, y = 1.045, label = "Low Neuroticism: b = -.12 [-.17, -.07], p < .001",
angle = -33.8)
## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.
A linear mixed model using restricted maximum likelihood was used to test whether the association of daily stress on daily mood is moderated by baseline neuroticism scores. All predictors were included as fixed effects and a random intercept by participant was included. Visual diagnostics showed that mood was not normally distributed, so it was square root transformed, referred to as transformed mood hereafter. Additionally, two outliers were identified on the residuals and these were excluded leaving a total of 181 observations from 50 people.
There was a significant daily stress x neuroticism interaction, such that higher neuroticism scores were associated with a more negative association between daily stress and same day transformed mood, b [95% CI] = -0.02 [-0.04, 0.00], p = .017. To help interpret the interaction, the simple slopes of daily stress on transformed mood were tested and graphed at low (M - 1 SD) and high (M + 1 SD) values of neuroticism (see Figure XX). The results revealed that there was a significant, negative association between daily stress and same day transformed mood at high, but not low levels of neuroticism.
Sensitivity analyses not excluding outliers revealed that results remained unchanged.
4 Continuous x Categorical Interactions in R
Continuous x Categorical interactions are conducted much as continuous x continuous interactions. Typically with continuous x categorical interactions, simple slopes for the continuous variable are calculated at all levels of the categorical variable.
Here we will work with a three level, dummy coded variable, stress with low being the reference and two dummy codes, one for mid and one for high stress. The outcome is transformed mood we used earlier.
Factor variables in interactions do not work currently with modelTest(), so if we wanted to use it, we’d need to manually dummy code the categorical variable. The results are identical.
## Parameters and CIs are based on REML,
## but modelTests requires ML not REML fit for comparisons,
## and these are used in effect sizes. Refitting.
## Term Est Type
## 1: (Intercept) 1.27*** [ 0.95, 1.59] Fixed Effects
## 2: neuroticism -0.01 [-0.05, 0.03] Fixed Effects
## 3: neuroticism:stressHigh -0.08 [-0.20, 0.03] Fixed Effects
## 4: neuroticism:stressMid -0.03 [-0.08, 0.03] Fixed Effects
## 5: stressHigh 0.38 [-0.56, 1.32] Fixed Effects
## 6: stressMid 0.12 [-0.30, 0.54] Fixed Effects
## 7: sd_(Intercept)|ID 0.19 Random Effects
## 8: sigma 0.29 Random Effects
## 9: Model DF 8 Overall Model
## 10: N (Groups) ID (50) Overall Model
## 11: N (Observations) 183 Overall Model
## 12: logLik -51.01 Overall Model
## 13: AIC 118.02 Overall Model
## 14: BIC 143.70 Overall Model
## 15: Marginal R2 0.11 Overall Model
## 16: Marginal F2 0.12 Overall Model
## 17: Conditional R2 0.37 Overall Model
## 18: Conditional F2 0.59 Overall Model
## 19: neuroticism (Fixed) 0.01/0.00, p = .638 Effect Sizes
## 20: stressMid (Fixed) 0.00/0.01, p = .578 Effect Sizes
## 21: stressHigh (Fixed) 0.01/0.00, p = .419 Effect Sizes
## 22: neuroticism:stressMid (Fixed) 0.01/0.01, p = .306 Effect Sizes
## 23: neuroticism:stressHigh (Fixed) 0.02/0.01, p = .144 Effect Sizes
## Term Est Type
4.1 Plotting
These interactions are not significant, so often we would simply drop the interaction term and report main effects. For the sake of example only, we will plot anyways.
With continuous x categorical interactions, the easiest approach is to plot the simple slope of the continuous variable by the categorical one as shown in the following.
These interactions are not significant, so often we would simply drop the interaction term and report main effects. For the sake of example only, we will continue to look at simple effects anyways.
When working with models that have interactions, a common aid to interpretation is to test the simple effects / slopes from the model. For example, previously we graphed the association between neuroticism and transformed mood at each level of the categorical stress3 variable. However, we cannot tell from the graph whether neuroticism is significantly associated with transformed mood for any specific level of stress. To answer that, we need to test the simple slope of neuroticism at specific values of neuroticism. Our default model does actually give us one simple slope: it is the simple slope of neuroticism when stress3 = low$, but we might want more.
We will use the emtrends() function from the emmeans package to test the simple slopes.
The emtrends() function take a model as its first argument, then the variable that you want to calculate a simple slope for, here stress, the argument at requires a list of specific values of the moderator, and then we tell it how we want degrees of freedom calculated (note this only applies to lmer models). We store the results in an R object, mem and then call summary() to get a summary table. The infer = TRUE argument is needed to summary() if you want p-values.
mem <- emtrends(m, var = "neuroticism",
at = list(stress3 = c("Low", "Mid", "High")),
lmer.df = "satterthwaite")
summary(mem, infer=TRUE)
The relevant parts of the output, for us, are the columsn for neuroticism.trend which are the simple slopes, the values of stress3 which tell us at what values of stress3 we have calculated simple slopes, the confidence intervals, lower.CL and upper.CL, 95% by default, and the p-value. From these results, we can see that neuroticism is not significantly associated with transformed mood for any stress category, although its closest for high stress.
5 Categorical x Categorical Interactions in R
Categorical x Categorical interactions are conducted comparably, although more contrasts / simple effect follow-ups are possible.
Here we will work with a three level, dummy coded variable, stress with low being the reference and two dummy codes, one for mid and one for high stress. The outcome is transformed mood we used earlier. We also work with a three level neuroticism variable.
With categorical x categorical interactions, visreg() produces OK but not great figures as shown in the following. We can see the means of transformed mood for all 9 cells (the cross of 3 level of neuroticism x 3 levels of stress).
When working with two categorical interactions (or, as an aside, with a categorical predictor with >2 levels where you want to test various group differences), the emmeans() function from the emmeans package is helpful. We can get the means of each neuroticism group by stress group and get confidence intervals and p-values. These p-values test whether each mean is different from zero, by default.
em <- emmeans(m, "neuro3", by = "stress3",
lmer.df = "satterthwaite")
summary(em, infer = TRUE)
A nice plot, with confidence intervals for the fixed effects, can be obtained by using the emmip() function from the emmeans package. It takes as input the results from emmeans(), not the lmer() model results directly. Here is a simple plot showing the categorical interactions. Note that with this approach, you could basically fit the same model(s) that you would with a repeated measures or mixed effects ANOVA model, with the advantage that LMMs do not require balanced designs and allow both categorical and continuous predictors (e.g., you could include continuous covariates easily). GLMs and (G)LMMs can do everything that t-tests and various ANOVAs can, but with greater flexibility.
If you want pairwise comparisons, you can get all possible pairwise comparisons between neuroticism levels by stress or stress levels by neuroticism using the pairs() function.
## pairwise comparisons of neuroticism by stress
pairs(em, by = "stress3")
## stress3 = Low:
## contrast estimate SE df t.ratio p.value
## Low - Mid -0.0310 0.117 89.9 -0.264 0.9625
## Low - High 0.0566 0.120 106.3 0.471 0.8851
## Mid - High 0.0875 0.141 105.8 0.622 0.8082
##
## stress3 = Mid:
## contrast estimate SE df t.ratio p.value
## Low - Mid -0.0501 0.107 75.3 -0.469 0.8861
## Low - High 0.1502 0.102 76.8 1.476 0.3079
## Mid - High 0.2004 0.114 70.9 1.751 0.1938
##
## stress3 = High:
## contrast estimate SE df t.ratio p.value
## Low - Mid -0.3432 0.217 173.9 -1.581 0.2568
## Low - High 0.1844 0.187 171.9 0.988 0.5853
## Mid - High 0.5276 0.199 170.6 2.657 0.0234
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 3 estimates
## pairwise comparisons of stress by neuroticism
pairs(em, by = "neuro3")
## neuro3 = Low:
## contrast estimate SE df t.ratio p.value
## Low - Mid 0.0746 0.0783 173 0.953 0.6073
## Low - High 0.3247 0.1493 160 2.175 0.0786
## Mid - High 0.2501 0.1480 157 1.690 0.2122
##
## neuro3 = Mid:
## contrast estimate SE df t.ratio p.value
## Low - Mid 0.0554 0.1070 174 0.518 0.8626
## Low - High 0.0125 0.1709 156 0.073 0.9971
## Mid - High -0.0429 0.1598 145 -0.269 0.9610
##
## neuro3 = High:
## contrast estimate SE df t.ratio p.value
## Low - Mid 0.1682 0.1056 170 1.594 0.2510
## Low - High 0.4526 0.1389 169 3.257 0.0039
## Mid - High 0.2843 0.1146 150 2.480 0.0377
##
## Degrees-of-freedom method: satterthwaite
## P value adjustment: tukey method for comparing a family of 3 estimates
You can also get custom contrasts. For example, if we wanted to compare high neuroticism to the average of low and mid neuroticism at each level of stress (H1) and secondly see if low and mid neuroticism differ (H2). The list gives the specific contrast weights, which are directly applied to the means we saw earlier from emmeans().
You could get even more specific hypotheses about group differences, if you wanted using by = NULL all the means are in a row of 9 and we can apply weights to each. For example, here we contrast the average of low neuroticism at low stress and low neuroticism at mid stress to high neuroticism and high stress (H1). For more examples of weighting schemes for many different possible specific contrasts, wee: https://stats.idre.ucla.edu/spss/faq/how-can-i-test-contrasts-and-interaction-contrasts-in-a-mixed-model/ .
